{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calibrate Pi Pulses\n",
    "*Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline\n",
    "This tutorial introduces how to calibrate a $\\pi$ pulse by varying the amplitude of the drive pulse. The outline of this tutorial is as follows:\n",
    "- Introduction\n",
    "- Preparation\n",
    "- Define the system Hamiltonian\n",
    "- Sweep amplitudes\n",
    "- Cosine regression\n",
    "- Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "Calibrating $\\pi$ pulses is one of the most fundamental operations in quantum computing, because one of the most fundamental gates, the X gate, requires a $\\pi$ pulse input onto the X channel. Further, it also serves an important role in calibrating actual hardware. Thus, this tutorial will demonstrate how to calibrate a $\\pi$ pulse using Quanlse."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "After you have successfully installed Quanlse, you could run the Quanlse program below following this tutorial. To run this particular tutorial, you would need to import the following packages from Quanlse and other commonly-used Python libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Hamiltonian module\n",
    "from Quanlse.Utils import Hamiltonian as qham \n",
    "\n",
    "# Import related packages\n",
    "from Quanlse.Utils.Operator import duff, driveX\n",
    "from Quanlse.Utils.Waveforms import gaussian\n",
    "\n",
    "# Import simulator interface for Quanlse Cloud Service\n",
    "from Quanlse.remoteSimulator import remoteSimulatorRunHamiltonian as runHamiltonian\n",
    "\n",
    "# Import numpy\n",
    "from numpy import linspace, pi, dot, array, cos\n",
    "\n",
    "# Import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Import curve_fit function from scipy\n",
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define the system Hamiltonian\n",
    "\n",
    "In the field of quantum control, it is a common practice to describe a quantum system with its Hamiltonian. Generally, a system Hamiltonian consists of two terms, the time-independent and the time-dependent terms:\n",
    "\n",
    "$$\n",
    "\\hat{H}_{\\rm total}(t) = \\hat{H}_{\\rm drift} + \\hat{H}_{\\rm ctrl }(t) .\n",
    "$$\n",
    "\n",
    "\n",
    "The users could easily define the Hamiltonian for a multi-qubit system using the `Hamiltonian` module in Quanlse. First, we will use the `Hamiltonian` module to initialize a system Hamiltonian. We start with a single-qubit system with three energy levels. The system Hamiltonian can be written as:\n",
    "\n",
    "$$\n",
    "\\hat{H} = \\alpha_q \\hat{a}^{\\dagger}\\hat{a}^{\\dagger}\\hat{a}\\hat{a} + \\frac{1}{2} c(t) \\cos(\\phi) (\\hat{a}+\\hat{a}^{\\dagger}).\n",
    "$$\n",
    "\n",
    "Here, the $\\alpha_q$ is the anharmonicity between the two lowest transisition energies. $c(t)$ indicates the pulse envelope function; and $\\phi$ is the pulse phase. $\\hat{a}^{\\dagger}=|1\\rangle\\langle 0|+\\sqrt{2}|2\\rangle\\langle 1|$ and $\\hat{a}=|0\\rangle\\langle 1|+\\sqrt{2}|1\\rangle\\langle 2|$ are respectively the creation and annihilation operators.\n",
    "\n",
    "Here, we will demonstrate how to define such a Hamiltonian using Quanlse. We will first initialize the Hamiltonian dictionary using the following code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ham = qham.createHam(title=\"example\", dt=0.2, qubitNum=1, sysLevel=3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above `createHam()` function returns an empty Hamiltonian dictionary. Its parameters include a user-defined title, sampling period, qubit number, and the system's energy levels to consider.\n",
    "\n",
    "Then we could add terms to the empty Hamiltonian dictionary using the two functions below. The function `addDrift()` adds drift operators to the Hamiltonian while the `addControl()` function adds the operators associated with the control pulses. Both functions require a `Hamiltonian` dictionary, a user-defined name, the qubit(s) index(es) which the term is acting upon, the according operators (we have conveniently provided the `Operator` module which includes many commonly-used operators), and the amplitude (only for the drift term):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "alphaq = - 0.22 * (2 * pi)  # unit is GHz\n",
    "qham.addDrift(ham, \"drift\", onQubits=0, matrices=duff(3), amp=alphaq)\n",
    "qham.addControl(ham, \"ctrl\", onQubits=0, matrices=driveX(3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, the user could use the `printHam()` function to display the properties of the Hamiltonian. The `printHam()` function takes a Hamiltonian dictionary:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qham.printHam(ham)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we could conveniently use `Operator`'s method `duff(n)` to define the $n$-dimensional $\\hat{a}^{\\dagger}\\hat{a}^{\\dagger}\\hat{a}\\hat{a}$, and `driveX(n)` to define the $n$-dimensional $\\frac{1}{2}(\\hat{a}+\\hat{a}^{\\dagger})$. After appending the control term to the Hamiltonian, we need to add the effective pulse:\n",
    "\n",
    "$$\n",
    "c(t) = A e^{-(t-\\tau)^2/2\\sigma^2}.\n",
    "$$\n",
    "\n",
    "We achieve this by using the `setWave()` function. The `setWave()` function takes a Hamiltonian dictionary, the name of the term in the Hamiltonian, waveform (Quanlse supports multiple waveforms' definitions), the parameters needed to define the wave, and lastly, the initial time and the duration of the wave."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qham.setWave(ham, \"ctrl\", f=\"gaussian\", para={\"a\": 1, \"tau\":10, \"sigma\":3}, t0=0, t=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we have just defined a complete quantum system and the parameters regarding controlling the system. We can visualize the pulse using the provided `plotWaves()` function. The `plotWaves()` function plots the pulses by taking a Hamiltonian dictionary and the according terms' names. The function also includes an optional bool parameter `dark`, which enables a dark-themed mode. Moreover, the user can use the `color` parameter to specify colors for individual pulses (the colors will repeat if there are more pulses than colors)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "qham.plotWaves(ham, \"ctrl\", dark=True, color=['mint'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can use the `simulate()` function to simulate the evolution, and obtain the unitary matrix of the system evolution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = qham.simulate(ham, recordEvolution=False)\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sweep amplitudes\n",
    "\n",
    "With fixed pulse duration $t_g$, we can sweep the pulse's amplitudes $a$, and find the amplitude $a_{\\pi}$ of the according $\\pi$ pulse.\n",
    "\n",
    "We first create a list of 200 points between -1 and 1, representing the pulse's amplitudes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initilize the pulse's amplitudes\n",
    "alist = linspace(-1.0, 1.0, 200)\n",
    "pop0_list = []\n",
    "pop1_list = []\n",
    "pop2_list = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we can obtain the according population for each state by simulating the evolution of the Hamiltonian defined in the previous section. The calculation usually takes a long time to process on local devices; however, we provide a cloud computing service that could speed up this process significantly. To use Quanlse Cloud Service, the users can get a token from http://quantum-hub.baidu.com and submit the job onto Quanlse's server. Note that Quanlse supports the submission of batches of job, which could further optimize the allocation of resources."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calibrate a Pi Pulse\n",
    "jobList = []\n",
    "for a in alist:\n",
    "    # Configure pulse parameters\n",
    "    jobTemp = []\n",
    "    jobTemp.append(qham.makeWaveData(ham, \"ctrl\", f=gaussian, para={\"a\": a, \"tau\": 10, \"sigma\": 3}, t0=0, t=20))\n",
    "    # Run similator\n",
    "    jobList.append(jobTemp)\n",
    "\n",
    "# Import Define class and set the token\n",
    "# Please visit http://quantum-hub.baidu.com\n",
    "from Quanlse import Define\n",
    "Define.hubToken = \"\"\n",
    "\n",
    "# Submit batch jobs to Quanlse Cloud Service\n",
    "resultList = runHamiltonian(ham, jobList=jobList)\n",
    "\n",
    "# Calculate populations\n",
    "for result in resultList:\n",
    "    final_state = dot(result[\"unitary\"], array([1, 0, 0], dtype=complex))\n",
    "    pop0_list.append(abs(final_state[0])**2)\n",
    "    pop1_list.append(abs(final_state[1])**2)\n",
    "    pop2_list.append(abs(final_state[2])**2)\n",
    "\n",
    "# Plot graph\n",
    "plt.plot(alist, pop0_list, label=\"Ground state\")\n",
    "plt.plot(alist, pop1_list, label=\"1st excited state\")\n",
    "plt.plot(alist, pop2_list, label=\"2nd excited state\")\n",
    "plt.xlabel(\"Amplitude\")\n",
    "plt.ylabel(\"Population of different states\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cosine regression\n",
    "\n",
    "Now, we have a series of discrete points; however, we need to fit those points with a cosine function in order to find the amplitude of the $\\pi$ pulse. To fit the resulting $|0\\rangle$ population, we use the `optimize.curve_fit()` method in `Scipy`. We first define the following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def fit_function(x_values, y_values, init_params):\n",
    "    def fit_func(x, A, B, period, phi):\n",
    "        return A * cos(2 * pi * x / period - phi) + B\n",
    "    fitparams, _ = curve_fit(fit_func, x_values, y_values, init_params, bounds=(0, [2.0, 2.0, 2.0, 2.0]))\n",
    "    y_fit = fit_func(x_values, *fitparams)\n",
    "    return fitparams, y_fit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we run the regression function to obtain the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fit_params, y_fit = fit_function(alist, pop0_list, [0.5, 0.5, 0.8, 0])\n",
    "\n",
    "# Plot graph\n",
    "plt.scatter(alist, pop0_list, label=\"Samples\")\n",
    "plt.plot(alist, y_fit, color=\"red\", label=\"Fit curve\")\n",
    "plt.xlabel(\"Amplitude\")\n",
    "plt.ylabel(\"Population of ground state\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "print(f\"Period is {fit_params[2]}\")\n",
    "print(f\"Pi pulse amplitude is {fit_params[2] / 2}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By the cosine regression, we have identified the corresponding amplitude of the $\\pi$ pulse."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "source": [
    "## Summary\n",
    "After reading this tutorial on calibrating $\\pi$ pulses, the users are encouraged to try parameter values different from this tutorial to obtain the optimal result."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
